In [2]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as dates
from mpl_toolkits.mplot3d import Axes3D
from datetime import datetime
import statsmodels.formula.api as sm
from scipy.optimize import fsolve
from scipy import stats
import random
In [3]:
dateparser=lambda x: pd.datetime.strptime(x,'%d/%m/%Y')
df = pd.read_csv('Data for Jupyter Lab.csv', sep=',',decimal='.', parse_dates=['Date_NORWAY','Date_SWEDEN','Date_LIBOR'], date_parser= dateparser)
In [4]:
df
Out[4]:
In [122]:
basic_statistics = pd.DataFrame()
basic_statistics['Norway']= [df.LN_NORWAY.mean(), df.LN_NORWAY.std(), stats.kurtosis(df.LN_NORWAY[1:-1])]
basic_statistics['Sweden']= [df.LN_SWEDEN.mean(), df.LN_SWEDEN.std(),stats.kurtosis(df.LN_SWEDEN[1:-1])]
basic_statistics['LIBOR']= [df.LN_LIBOR.mean(), df.LN_LIBOR.std(), stats.kurtosis(df.LN_LIBOR[1:-1])]
basic_statistics.rename(index = {0:"Mean", 1:"Standard Deviation", 2:"Kurtosis"}, inplace = True)
basic_statistics
Out[122]:
In [59]:
plt.grid(b=True)
plt.xlabel('Date')
plt.ylabel('Price (USD/NOK)')
plt.title("Price path for\nNorwegian exchange rate")
plt.fill_between(list(df.Date_NORWAY), df.Price_NORWAY, color="red", alpha=0.2)
Out[59]:
In [123]:
plt.grid(b=True)
plt.xlabel('Date')
plt.ylabel('Price (SEK)')
plt.title("Price path for\nSwedish stock market")
plt.fill_between(list(df.Date_SWEDEN), df.Price_SWEDEN, color="gold", alpha=0.2)
Out[123]:
In [6]:
plt.grid(b=True)
plt.xlabel('Date')
plt.ylabel('Annual Yield (%)')
plt.title("Yield path for\nLIBOR 12-month rate in British pounds")
plt.fill_between(list(df.Date_LIBOR), df.Price_LIBOR, color="blue", alpha=0.2)
Out[6]:
In [7]:
plt.figure(figsize=(24,5))
plt.grid(b=True)
plt.xlabel('Date')
plt.ylabel('Price change\n(%)')
plt.title("Natural-logs for\nNorwegian exchange rate")
plt.plot(df.Date_NORWAY, df.LN_NORWAY, color="red",linewidth=0.5)
plt.gca().set_yticklabels(['{:.0f}'.format(x*100) for x in plt.gca().get_yticks()])
plt.show()
In [8]:
plt.figure(figsize=(24,5))
plt.grid(b=True)
plt.xlabel('Date')
plt.ylabel('Price change\n(%)')
plt.title("Natural-logs for\nSwedish stock market")
plt.plot(df.Date_SWEDEN, df.LN_SWEDEN, color="gold",linewidth=0.5)
plt.gca().set_yticklabels(['{:.0f}'.format(x*100) for x in plt.gca().get_yticks()])
plt.show()
In [9]:
plt.figure(figsize=(24,5))
plt.grid(b=True)
plt.xlabel('Date')
plt.ylabel('Yield change\n(%)')
plt.title("Natural-logs for\nLIBOR 12-month rate in British pounds")
plt.plot(df.Date_LIBOR, df.LN_LIBOR, color="blue",linewidth=0.5)
plt.gca().set_yticklabels(['{:.0f}'.format(x*100) for x in plt.gca().get_yticks()])
plt.show()
In [10]:
print("Here we consider if the daily price-changes for our three markets are Normally distributed.\nTo do this we use two statistical tests.")
print("\nHere we show the resulting p-values.\n\nAs we can see, the odds of these markets returns being Normally distributed are practically zero.\n")
import warnings
warnings.simplefilter("ignore")
Statistical_Tests=[[stats.kstest(list(df.LN_NORWAY[1:-1]), 'norm').pvalue,
stats.kstest(list(df.LN_SWEDEN[1:-1]), 'norm').pvalue,
stats.kstest(list(df.LN_LIBOR[1:-1]), 'norm').pvalue] ,
[stats.shapiro(list(df.LN_NORWAY[1:-1]))[1],
stats.shapiro(list(df.LN_SWEDEN[1:-1]))[1],
stats.shapiro(list(df.LN_LIBOR[1:-1]))[1]]]
Statistical_Tests = pd.DataFrame(Statistical_Tests, columns = ['Norway', 'Sweden', 'LIBOR'], index = ['Kolmogorov-Smirnov', 'Shapiro-Wilk'])
Statistical_Tests
Out[10]:
In [10]:
print("To begin the multifractal analysis, we must first define our Xt variables.\nIt is useful to have a function for this.")
Xt= lambda P: np.log(P) - np.log(P[0])
df['Xt_NORWAY']= Xt(df.Price_NORWAY)
df['Xt_SWEDEN']= Xt(df.Price_SWEDEN)
df['Xt_LIBOR']= Xt(df.Price_LIBOR)
In [11]:
df[['Price_NORWAY','LN_NORWAY','Xt_NORWAY','Price_SWEDEN','LN_SWEDEN','Xt_SWEDEN','Price_LIBOR','LN_LIBOR','Xt_LIBOR']]
Out[11]:
In [13]:
plt.grid(b=True)
plt.xlabel('t\n(time in days)')
plt.ylabel('Xt \n(natural-log change since time 0)')
plt.title("Growth over time for Norwegian exchange rate")
plt.fill_between(list(df.t), df.Xt_NORWAY, color="red", alpha=0.2)
Out[13]:
In [14]:
plt.grid(b=True)
plt.xlabel('t\n(time in days)')
plt.ylabel('Xt \n(natural-log change since time 0)')
plt.title("Growth over time for Swedish stock market")
plt.fill_between(list(df.t), df.Xt_SWEDEN, color="gold", alpha=0.2)
Out[14]:
In [15]:
plt.grid(b=True)
plt.xlabel('t\n(time in days)')
plt.ylabel('Xt \n(natural-log change since time 0)')
plt.title("Growth over time for British LIBOR")
plt.fill_between(list(df.t), df.Xt_LIBOR, color="blue", alpha=0.2)
Out[15]:
In [12]:
print("We will be using a HIGHLY COMPOSITE NUMBER - 7560 - as our number of observations. \nThis is because this number has a lot of factors (numbers that can divide 7560 and yield a whole number).\n \nWe will call our time increments 'delta_t'.\n \nThe number of time increments - delta_t's - that we use in the analysis is the same as the number of factors. ")
delta_t=(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 18, 20, 21, 24, 27, 28, 30, 35, 36, 40, 42, 45, 54, 56, 60, 63, 70, 72, 84, 90, 105, 108, 120, 126, 135, 140, 168, 180, 189, 210, 216, 252, 270, 280, 315, 360, 378, 420, 504, 540, 630, 756, 840, 945, 1080, 1260, 1512, 1890, 2520, 3780, 7560)
In [13]:
print("The number of delta_t's is: " + str(len(delta_t)))
In [14]:
print("Here are all the delta_t's (the factors of 7560):\n")
for i in range (0,int(len(delta_t)/2)):
print("Factor " + str(i+1) + ": 7560 divided by " + str(delta_t[i]) + " is " + str(int(7560/delta_t[i])))
print("\n...\n \netc. \n \nAfter factor 32, the next factors are simply the results of the \ndivisions (90, 105, 108 etc. up to 7560), for a total of 64 factors.")
In [15]:
print("Let us now choose the values of q - the statistical moments - that we will use for our partition function.")
print("\n(Recall that q=1 is the mean, q=2 is the variance, q=4 is the kurtosis etc.")
print("Well... almost... in fact, only q=1 is as the mean is correct. The other 'raw' moments need to be normalized.\nBut we don't need that distinction here.)")
print("\nMandelbrot et al (1997) chose approximately 99 values for q, mostly concentrated around the low values between 0 and 5.\nWe shall try to do something similar.")
q=[0.01, 0.1, 0.2, 0.3, 0.4,
0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
1.0, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3, 1.35, 1.4, 1.45,
1.5, 1.55, 1.6, 1.65, 1.7, 1.75,
1.8, 1.81, 1.82, 1.83, 1.84, 1.85, 1.86, 1.87, 1.88, 1.89,
1.9, 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97, 1.98, 1.985, 1.99,
1.991, 1.992, 1.993, 1.994, 1.995, 1.996, 1.997, 1.998, 1.999,
2.0,
2.001, 2.002, 2.003, 2.004, 2.005, 2.006, 2.007, 2.008, 2.009,
2.01, 2.015, 2.02, 2.025,
2.03, 2.04, 2.05, 2.06, 2.07, 2.08, 2.09,
2.1, 2.15, 2.2, 2.25, 2.3, 2.35, 2.4, 2.45, 2.5,
2.6, 2.7, 2.8, 2.9,
3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9,
4.0, 4.5,
5.0, 6.0, 7.0, 8.0, 9.0, 10.0,
12.5, 15.0, 17.5, 20.0, 22.5, 25.0, 27.5,
30.0]
print("\nThe number of q's that we will use is " + str(len(q)) + ".")
In [16]:
plt.figure(figsize=(12,3))
plt.grid(b=True)
plt.title("The distribution of q's used for the Partition Function (Sq(t,delta_t))")
plt.xlabel("Index of q\n(Total " + str(len(q)) + ")")
plt.ylabel("q\n(the 'raw' moment)")
plt.plot(q, marker='x')
Out[16]:
In [17]:
print("q is:\n\n" + str(q))
In [18]:
print("The number of q's between 1 and 3 is: " + str(sum(float(num) <= 3 and float(num) >=1 for num in q)) +"\n")
print("The number of q's less than 1 is: " + str(sum(float(num) <= 1 for num in q)))
print("The number of q's greater than 3 is: " + str(sum(float(num) >= 3 for num in q)) +"\n")
print("Total q's: " + str(len(q)))
In [19]:
print("The next step is to figure out the partition function - Sq(T,delta_t) - for different values of delta_t and q.")
print("The partition function Sq(T,delta_t) will be used to estimate the scaling function tau(q).")
def partition_function(SIGMA, DELTA, XT, Q):
print("Calculating the partition function...\nThis step will take quite a while... so strap yourself in...\n")
SIGMA=[[0 for x in range(len(DELTA))] for y in range(len(Q))]
for k in range (0, len(Q)):
if k%30==0:
print("calculating i=" + str(k) + ' out of ' + str(len(Q)-1))
for j in range (0,len(DELTA)):
for i in range (0,len(XT)-1):
if i < int(len(XT)/DELTA[j]):
SIGMA[k][j]=SIGMA[k][j] + abs(XT[i*DELTA[j]+DELTA[j]]-XT[i*DELTA[j]])**Q[k]
SIGMA=pd.DataFrame(SIGMA)
for i in range (0,len(Q)):
SIGMA.rename(index={SIGMA.index[i]:Q[i]}, inplace=True)
for i in range (len(DELTA)-1,-1,-1):
SIGMA.rename(columns={SIGMA.columns[i]:DELTA[i]}, inplace=True)
print("Done! Your partition function is ready!\n")
return SIGMA
print("\nHere we have defined our partition function in the form: \npartition_function(SIGMA, DELTA, XT, Q).")
In [20]:
print("Here, we will now begin calculating the partition functions for all our markets.\n")
partition_NORWAY=[[]]
partition_SWEDEN=[[]]
partition_LIBOR=[[]]
partition_NORWAY=partition_function(partition_NORWAY, delta_t, df.Xt_NORWAY, q)
partition_SWEDEN=partition_function(partition_SWEDEN, delta_t, df.Xt_SWEDEN, q)
partition_LIBOR=partition_function(partition_LIBOR, delta_t, df.Xt_LIBOR, q)
print("\nDone! You can now proceed to estimating the scaling function - tau(q).")
In [21]:
for i in range (0, len(q)):
plt.plot(np.log(delta_t), np.log(list(partition_NORWAY.iloc[i])/partition_NORWAY[1][q[i]]), color="red",linewidth=0.5, label=str(q[i]))
plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
ncol=6, mode="expand", borderaxespad=0.)
plt.xlabel('ln (delta_t)\n(the natural log of time increments)')
plt.ylabel('ln ( Sq(delta_t) )\n(the natural log of the partition function)')
plt.show()
In [22]:
for i in range (0, len(q)):
plt.plot(np.log(delta_t), np.log(list(partition_SWEDEN.iloc[i])/partition_SWEDEN[1][q[i]]), color="gold",linewidth=0.5, label=str(q[i]))
plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
ncol=6, mode="expand", borderaxespad=0.)
plt.xlabel('ln (delta_t)\n(the natural log of time increments)')
plt.ylabel('ln ( Sq(delta_t) )\n(the natural log of the partition function)')
plt.show()
In [23]:
for i in range (0, len(q)):
plt.plot(np.log(delta_t), np.log(list(partition_LIBOR.iloc[i])/partition_LIBOR[1][q[i]]), color='blue',linewidth=0.5, label=str(q[i]))
plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
ncol=6, mode="expand", borderaxespad=0.)
plt.xlabel('ln (delta_t)\n(the natural log of time increments)')
plt.ylabel('ln ( Sq(delta_t) )\n(the natural log of the partition function)')
plt.show()
In [ ]:
In [24]:
for i in range (0, len(q)):
if q[i]%0.25==0 and q[i]<6 and q[i]>1 and q[i]!=3.5 and q[i]!=4.5:
plt.plot(np.log(delta_t), np.log(list(partition_NORWAY.iloc[i])/partition_NORWAY[1][q[i]]), color="red", linewidth=1.5, label='q='+str(q[i]))
plt.legend(bbox_to_anchor=(0., -0.45, 1., .102), loc=3,
ncol=5, mode="expand", borderaxespad=0., title="Raw moments")
plt.title("Parition function for Norway")
plt.xlabel('ln (delta_t)\n(the natural log of time increments)')
plt.ylabel('ln ( Sq(delta_t) )\n(the natural log of the partition function)')
plt.show()
In [25]:
for i in range (0, len(q)):
if q[i]%0.25==0 and q[i]<6 and q[i]>1 and q[i]!=3.5 and q[i]!=4.5:
plt.plot(np.log(delta_t), np.log(list(partition_SWEDEN.iloc[i])/partition_SWEDEN[1][q[i]]), color="gold",linewidth=1.5, label='q='+str(q[i]))
plt.legend(bbox_to_anchor=(0., -0.45, 1., .102), loc=3,
ncol=5, mode="expand", borderaxespad=0., title="Raw moments")
plt.title("Parition function for Sweden")
plt.xlabel('ln (delta_t)\n(the natural log of time increments)')
plt.ylabel('ln ( Sq(delta_t) )\n(the natural log of the partition function)')
plt.show()
In [26]:
for i in range (0, len(q)):
if q[i]%0.25==0 and q[i]<6 and q[i]>1 and q[i]!=3.5 and q[i]!=4.5:
plt.plot(np.log(delta_t), np.log(list(partition_LIBOR.iloc[i])/partition_LIBOR[1][q[i]]), color="blue",linewidth=1.5, label='q='+str(q[i]))
plt.legend(bbox_to_anchor=(0., -0.45, 1., .102), loc=3,
ncol=5, mode="expand", borderaxespad=0., title="Raw moments")
plt.title("Parition function for LIBOR")
plt.xlabel('ln (delta_t)\n(the natural log of time increments)')
plt.ylabel('ln ( Sq(delta_t) )\n(the natural log of the partition function)')
plt.show()
In [27]:
tau_regression=pd.DataFrame(np.log(delta_t))
tau_regression['LN_DELTA']=pd.DataFrame(np.log(delta_t))
tau_regression['LN_T']=[np.log(7560) for x in range(len(delta_t))]
def scaling_function(TAU_Q, PARTITION, Q):
TAU_Q=[0 for x in range(len(Q))]
for i in range(0,len(Q)):
TAU_Q[i]=((sm.OLS(endog=np.log(list(PARTITION.iloc[i]/PARTITION[1][Q[i]])), exog = tau_regression[['LN_DELTA','LN_T']], missing='drop')).fit()).params[0]
TAU_Q=pd.DataFrame(TAU_Q)
return TAU_Q
print("Here we have defined how to find the estimated scaling function tau_q. \n\nWe are now ready to estimate tau_q for each market.")
In [28]:
tau_q_NORWAY=0
tau_q_SWEDEN=0
tau_q_LIBOR=0
tau_q_NORWAY = scaling_function(tau_q_NORWAY, partition_NORWAY, q)
tau_q_SWEDEN = scaling_function(tau_q_SWEDEN, partition_SWEDEN, q)
tau_q_LIBOR = scaling_function(tau_q_LIBOR, partition_LIBOR, q)
hypothetical_tau_q = [0.5*x -1 for x in q]
In [29]:
plt.grid(b=True)
plt.title("The scaling function for tau_q_NORWAY")
plt.xlabel("q\n(the 'raw' moment)")
plt.ylabel("tau(q)\n(the 'scaling function')")
axes = plt.gca()
axes.set_xlim([-0.5,8.5])
axes.set_ylim([-1,2.5])
plt.plot(q[0:111], tau_q_NORWAY[0:111], color="red",linewidth=3, marker="x", label="observed tau(q)")
plt.plot(q[0:111], hypothetical_tau_q[0:111], color="black", linestyle='--', label="hypothetical monofractal")
plt.legend()
plt.show()
In [30]:
plt.grid(b=True)
plt.title("The scaling function for tau_q_SWEDEN")
plt.xlabel("q\n(the 'raw' moment)")
plt.ylabel("tau(q)\n(the 'scaling function')")
axes = plt.gca()
axes.set_xlim([-0.5,8.5])
axes.set_ylim([-1,2.5])
plt.plot(q[0:111], tau_q_SWEDEN[0:111], color="gold",linewidth=3, marker="x", label="observed tau(q)")
plt.plot(q[0:111], hypothetical_tau_q[0:111], color="black", linestyle='--', label="hypothetical monofractal")
plt.legend()
plt.show()
In [31]:
plt.grid(b=True)
plt.title("The scaling function for tau_q_LIBOR")
plt.xlabel("q\n(the 'raw' moment)")
plt.ylabel("tau(q)\n(the 'scaling function')")
axes = plt.gca()
axes.set_xlim([-0.5,8.5])
axes.set_ylim([-1,2.5])
plt.plot(q[0:111], tau_q_LIBOR[0:111], color="blue",linewidth=3, marker="x", label="observed tau(q)")
plt.plot(q[0:111], hypothetical_tau_q[0:111], color="black", linestyle='--', label="hypothetical monofractal")
plt.legend()
plt.show()
In [32]:
print("Our intercept is somewhere between:")
print(str(tau_q_NORWAY[0][85])+ " for q=" +str(q[85]))
print(str(tau_q_NORWAY[0][86])+ " for q=" +str(q[86]))
print("\nTherefore NORWAY's H is " + str(round(1/q[86], 3)) + " < H < " + str(round(1/q[85], 3)))
In [33]:
print("Our intercept is somewhere between:")
print(str(tau_q_SWEDEN[0][35])+ " for q=" +str(q[35]))
print(str(tau_q_SWEDEN[0][36])+ " for q=" +str(q[36]))
print("\nTherefore Sweden's H is " + str(round(1/q[36], 3)) + " < H < " + str(round(1/q[35], 3)))
In [34]:
print("Our intercept is somewhere between:")
print(str(tau_q_LIBOR[0][26])+ " for q=" +str(q[26]))
print(str(tau_q_LIBOR[0][27])+ " for q=" +str(q[27]))
print("\nTherefore the UK LIBOR's H is " + str(round(1/q[27], 3)) + " < H < " + str(round(1/q[26], 3)))
In [35]:
max_q=105
print("Here we have defined the highest moment q which we will use for estimating the multifractal spectrum.")
In [36]:
print("We have seen that the partition function Sq(T,delta_t) for higher moments is less robust.\nSo in estimating our scaling function, we will restrict ourselves to lower moments.")
print("\nHere, we will only use moments up to q=" + str(q[max_q]) + " to estimate the scaling function tau(q) by non-linear regression.")
print("\nNext, using the parameters that we obtained from regression, we will perform\na Legendre transformation on tau(q) to estimate the multifractal spectrum f(a).")
def estimate_multifractal_spectrum(TAU_Q, Q, MIN_Q, MAX_Q):
TAU_Q_ESTIMATED = np.polyfit(Q[MIN_Q:MAX_Q], TAU_Q[MIN_Q:MAX_Q], 2)
F_A = [0 for x in range(len(q)-10)]
p = [0 for x in range(len(q)-10)]
a = TAU_Q_ESTIMATED[0][0]
b = TAU_Q_ESTIMATED[1][0]
c = TAU_Q_ESTIMATED[2][0]
for i in range(0, len(q)-10):
p[i] = 2*a*Q[i]+b
F_A[i] = ((p[i]-b)/(2*a))*p[i] - (a*((p[i]-b)/(2*a))**2 + b*((p[i]-b)/(2*a)) + c)
F_A = pd.DataFrame(F_A)
F_A.rename(columns={F_A.columns[0]:"f(a)"}, inplace=True)
F_A['p'] = p
print("Using the range of q's from " + str(Q[MIN_Q]) + " to " + str(Q[MAX_Q]) + ":")
print("The estimated parameters for tau(q) are: \n" + str(TAU_Q_ESTIMATED))
print("\nThus, the estimated parameters for f(a) are: \n" + str(1/(4*a)) + ", \n" + str((-2*b)/(4*a)) + ", \n"+ str((-4*a*c+b**2)/(4*a)))
return F_A
In [37]:
f_a_NORWAY = estimate_multifractal_spectrum(tau_q_NORWAY, q, 0, max_q)
print ("\nf_a_NORWAY[f(a)][0] is: " + str(f_a_NORWAY['f(a)'][0]))
In [38]:
f_a_SWEDEN = estimate_multifractal_spectrum(tau_q_SWEDEN, q, 0, max_q)
print ("\nf_a_SWEDEN[f(a)][0] is: " + str(f_a_SWEDEN['f(a)'][0]))
In [39]:
f_a_LIBOR = estimate_multifractal_spectrum(tau_q_LIBOR, q, 0, max_q)
print ("\nf_a_LIBOR[f(a)][0] is: " + str(f_a_LIBOR['f(a)'][0]))
In [40]:
plt.grid(b=True)
plt.title("The multifractal spectrum for f_a_NORWAY")
plt.xlabel("a\n(the Hölder exponent)")
plt.ylabel('f(a)\n(the multifractal spectrum)')
# axes = plt.gca()
# axes.set_xlim([0.2,0.7])
# axes.set_ylim([0,1.1])
plt.plot(f_a_NORWAY['p'],f_a_NORWAY['f(a)'], color="red",linewidth=5)
Out[40]:
In [41]:
plt.grid(b=True)
plt.title("The multifractal spectrum for f_a_SWEDEN")
plt.xlabel("a\n(the Hölder exponent)")
plt.ylabel('f(a)\n(the multifractal spectrum)')
# axes = plt.gca()
# axes.set_xlim([0.2,0.8])
# axes.set_ylim([0,1.1])
plt.plot(f_a_SWEDEN['p'],f_a_SWEDEN['f(a)'], color="gold",linewidth=5)
Out[41]:
In [42]:
plt.grid(b=True)
plt.title("The multifractal spectrum for f_a_LIBOR")
plt.xlabel("a\n(the Hölder exponent)")
plt.ylabel('f(a)\n(the multifractal spectrum)')
# axes = plt.gca()
# axes.set_xlim([0.2,0.8])
# axes.set_ylim([0,1.1])
plt.plot(f_a_LIBOR['p'],f_a_LIBOR['f(a)'], color="blue",linewidth=5)
Out[42]:
In [43]:
print("Using the parameters for our tau(q)'s that we estimated by non-linear regression, \nwe can calculate precise values for H:")
In [44]:
def estimate_H(TAU_Q, Q, MIN_Q, MAX_Q):
TAU_Q_ESTIMATED = np.polyfit(Q[MIN_Q:MAX_Q], TAU_Q[MIN_Q:MAX_Q], 2)
def f(x):
return TAU_Q_ESTIMATED[0][0]*x**2 + TAU_Q_ESTIMATED[1][0]*x + TAU_Q_ESTIMATED[2][0]
temp = fsolve(f, [0,4])
H = 1/temp
return H[0]
print("Here we have defined our function for precisely estimating H.")
In [45]:
H_NORWAY = estimate_H(tau_q_NORWAY, q, 0, max_q)
H_SWEDEN = estimate_H(tau_q_SWEDEN, q, 0, max_q)
H_LIBOR = estimate_H(tau_q_LIBOR, q, 0, max_q)
print("Here, we calculate our H-values to be:\n")
print("For NORWAY: H = " + str(H_NORWAY))
print("For Sweden: H = " + str(H_SWEDEN))
print("For LIBOR: H = " + str(H_LIBOR))
In [46]:
print("Thus, using our estimated multifractal spectra, we get that:")
print("\nNORWAY's estimated alpha zero = " + str(f_a_NORWAY['p'][0]))
print("\nSweden's estimated alpha zero = " + str(f_a_SWEDEN['p'][0]))
print("\nLIBOR's estimated alpha zero = " + str(f_a_LIBOR['p'][0]))
In [47]:
print("From this, it is easy to estimate the means and variances for the log-normal distribution:")
print("\nUsing: \n \n𝜆 = a / H \n \n and \n \n𝜎^2 = (2(𝜆-1))/ln[b]")
def lambda_mean(H,ALPHA):
LAMBDA = ALPHA/H
return LAMBDA
def sigma_variance(LAMBDA, B):
SIGMA_VARIANCE = (2*(LAMBDA-1))/np.log(B)
return SIGMA_VARIANCE
In [48]:
print("Assuming we will partition our multifractal cascade in two at each step, meaning b=2:")
b=2
lambda_NORWAY = lambda_mean(H_NORWAY, f_a_NORWAY['p'][0])
lambda_SWEDEN = lambda_mean(H_SWEDEN, f_a_SWEDEN['p'][0])
lambda_LIBOR = lambda_mean(H_LIBOR, f_a_LIBOR['p'][0])
sigma_NORWAY = sigma_variance(lambda_NORWAY, b)
sigma_SWEDEN = sigma_variance(lambda_SWEDEN, b)
sigma_LIBOR = sigma_variance(lambda_LIBOR, b)
print("\nFor NORWAY, we estimate\n 𝜆 = " + str(lambda_NORWAY) + "\n 𝜎^2 = " + str(sigma_NORWAY))
print("\nFor Sweden, we estimate\n 𝜆 = " + str(lambda_SWEDEN) + "\n 𝜎^2 = " + str(sigma_SWEDEN))
print("\nFor LIBOR, we estimate\n 𝜆 = " + str(lambda_LIBOR) + "\n 𝜎^2 = " + str(sigma_LIBOR))
In [49]:
print("Here, we can print our final estimated results.\n")
print("We are now ready to begin constructing our MMAR simulations!\n")
RESULTS = pd.DataFrame([0 for x in range(0,3)])
# RESULTS = pd.DataFrame([[H_NORWAY,H_SWEDEN,H_LIBOR],
# [f_a_NORWAY['p'][0],f_a_SWEDEN['p'][0], f_a_LIBOR['p'][0]],
# [lambda_NORWAY, lambda_SWEDEN, lambda_LIBOR]
# [sigma_NORWAY, sigma_SWEDEN, sigma_LIBOR]])
RESULTS['H'] = pd.DataFrame([H_NORWAY,H_SWEDEN,H_LIBOR])
RESULTS['alpha zero'] = pd.DataFrame([f_a_NORWAY['p'][0],f_a_SWEDEN['p'][0], f_a_LIBOR['p'][0]])
RESULTS['lambda'] = pd.DataFrame([lambda_NORWAY, lambda_SWEDEN, lambda_LIBOR])
RESULTS['sigma^2'] = pd.DataFrame([sigma_NORWAY, sigma_SWEDEN, sigma_LIBOR])
RESULTS.rename(index={RESULTS.index[0]:'Norway'}, inplace=True)
RESULTS.rename(index={RESULTS.index[1]:'Sweden'}, inplace=True)
RESULTS.rename(index={RESULTS.index[2]:'UK LIBOR'}, inplace=True)
RESULTS.drop(labels=0, axis=1, inplace=True)
RESULTS
RESULTS.round(decimals=3)
Out[49]:
In [50]:
RESULTS
Out[50]:
In [127]:
tau_q_ESTIMATED_NORWAY=0
tau_q_ESTIMATED_NORWAY=np.polyfit(q[0:max_q],tau_q_NORWAY[0:max_q],2)
for i in range(0,max_q):
plt.scatter(q[i], (tau_q_ESTIMATED_NORWAY[0][0]*(q[i]**2) + tau_q_ESTIMATED_NORWAY[1][0]*q[i] + tau_q_ESTIMATED_NORWAY[2][0]), s = 20, marker='x')
plt.plot(q[0:max_q], tau_q_NORWAY[0:max_q], color="red",linewidth=5, alpha = 0.5)
#plt.plot(q,tau_q_ESTIMATED_NORWAY, linewidth=1)
plt.grid(b=True)
plt.title("Estimated tau(q) for Norway")
plt.xlabel("q\n(the 'raw' moment)")
plt.axvline(1/H_NORWAY, color="black", linestyle='--', linewidth=0.9, label = "H = 1 / q\n = " + str(round(H_NORWAY, 3)))
plt.axhline(0, color="black", linestyle='--', linewidth=0.9)
plt.ylabel("tau(q)\n(the 'scaling function')")
plt.legend()
plt.show()
print("\nFor NORWAY, we estimate\n 𝜆 = " + str(lambda_NORWAY) + "\n 𝜎^2 = " + str(sigma_NORWAY))
In [128]:
tau_q_ESTIMATED_SWEDEN=0
tau_q_ESTIMATED_SWEDEN=np.polyfit(q[0:max_q],tau_q_SWEDEN[0:max_q],2)
for i in range(0,max_q):
plt.scatter(q[i], (tau_q_ESTIMATED_SWEDEN[0][0]*(q[i]**2) + tau_q_ESTIMATED_SWEDEN[1][0]*q[i] + tau_q_ESTIMATED_SWEDEN[2][0]), s = 20, marker='x')
plt.plot(q[0:max_q], tau_q_SWEDEN[0:max_q], color="gold",linewidth=5, alpha=0.5)
#plt.plot(q,tau_q_ESTIMATED_NORWAY, linewidth=1)
plt.grid(b=True)
plt.title("Estimated tau(q) for Sweden")
plt.xlabel("q\n(the 'raw' moment)")
plt.axvline(1/H_SWEDEN, color="black", linestyle='--', linewidth=0.9, label = "H = 1 / q\n = " + str(round(H_SWEDEN, 3)))
plt.axhline(0, color="black", linestyle='--', linewidth=0.9)
plt.ylabel("tau(q)\n(the 'scaling function')")
plt.legend()
plt.show()
print("\nFor SWEDEN, we estimate\n 𝜆 = " + str(lambda_SWEDEN) + "\n 𝜎^2 = " + str(sigma_SWEDEN))
In [129]:
tau_q_ESTIMATED_LIBOR=0
tau_q_ESTIMATED_LIBOR=np.polyfit(q[0:max_q],tau_q_LIBOR[0:max_q],2)
for i in range(0,max_q):
plt.scatter(q[i], (tau_q_ESTIMATED_LIBOR[0][0]*(q[i]**2) + tau_q_ESTIMATED_LIBOR[1][0]*q[i] + tau_q_ESTIMATED_LIBOR[2][0]), s = 20, marker='x')
plt.plot(q[0:max_q], tau_q_LIBOR[0:max_q], color="blue",linewidth=5, alpha=0.5)
#plt.plot(q,tau_q_ESTIMATED_NORWAY, linewidth=1)
plt.grid(b=True)
plt.title("Estimated tau(q) for LIBOR")
plt.xlabel("q\n(the 'raw' moment)")
plt.axvline(1/H_LIBOR, color="black", linestyle='--', linewidth=0.9, label = "H = 1 / q\n = " + str(round(H_LIBOR, 3)))
plt.axhline(0, color="black", linestyle='--', linewidth=0.9)
plt.ylabel("tau(q)\n(the 'scaling function')")
plt.legend()
plt.show()
print("\nFor LIBOR, we estimate\n 𝜆 = " + str(lambda_LIBOR) + "\n 𝜎^2 = " + str(sigma_LIBOR))
In [ ]:
In [55]:
pip install fbm
In [5]:
from fbm import FBM
from fbm import fbm, fgn, times
from fbm import MBM
from fbm import mbm, mgn, times
import math
In [127]:
def lognormal_cascade(k, v,ln_lambda, ln_theta):
k = k - 1
m0 = np.random.lognormal(ln_lambda,ln_theta)
m1 = np.random.lognormal(ln_lambda,ln_theta)
M = [m0, m1]
if (k >= 0):
d=[0 for x in range(0,2)]
for i in range(0,2):
d[i] = lognormal_cascade(k, (M[i]*v), ln_lambda, ln_theta)
v = d
return v
In [129]:
def MMAR(K, simulated_H, simulated_lambda, simulated_sigma, original_price_history, magnitude_parameter, GRAPHS):
# --- VARIABLES ---
# K
# adjust K depending on how many days you want to simulate (e.g. if K=13, you'll simulate 2^13=8192 days)
# simulated_H
# the Hurst parameter for the fBm process
# simulated_lambda
# the mean for the lognormal cascade
# simulated_sigma
# the variance for the lognormal cascade
# original_price_history
# the price history of the market under study (used for starting the prices from the right time!)
# magnitude_parameter
# adjust this up or down to change the range of price changes (e.g. if prices are swinging too wildly every day, then adjust this downwards)
# GRAPHS
# Boolean - either True or False - use True if you want the MMAR to simulate graphs for you
# --- MESSAGE ---
if GRAPHS == True:
print("Performing an MMAR simulation with parameters:\n\nH = " + str(simulated_H) + "\nlambda = " + str(simulated_lambda) + "\nsigma = " + str(simulated_sigma) + "\nfBm magnitude = " + str(magnitude_parameter)+ "\n")
# --- CASCADE ---
new_cascade = list(np.array(lognormal_cascade(k=K, v=1, ln_lambda = simulated_lambda, ln_theta = simulated_sigma)).flat)
if GRAPHS == True:
# plt.figure(figsize=(24,2))
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.title("Binomial Cascade")
plt.xlabel("Conventional time\n(days)")
plt.ylabel('"Mass"')
plt.plot(new_cascade, color="crimson", linewidth=0.5)
plt.show()
# --- TRADING TIME ---
tradingtime = 2**K*np.cumsum(new_cascade)/sum(new_cascade)
if GRAPHS == True:
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.yticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.title("Trading time\n(normalized)")
plt.xlabel("Conventional time\n(days)")
plt.ylabel('"Trading time"\n(\u03B8 (t), normalized)')
plt.plot(tradingtime, color="orangered")
# --- FBM (Fractional Brownian Motion) ---
new_fbm_class = FBM(n = 10*2**K+1, hurst = simulated_H, length = magnitude_parameter, method='daviesharte')
new_fbm_simulation = new_fbm_class.fbm()
if GRAPHS == True:
plt.figure(figsize=(24,2))
plt.xticks(np.arange(0, 10*2**(K)+1, 10*2**(K-3)))
plt.title("Fractional Brownian Motion")
plt.xlabel("t")
plt.ylabel('fBm (t)')
plt.plot(new_fbm_simulation, color="orange")
plt.show()
# --- MMAR XT's ---
simulated_xt_array = [0 for x in range(0, len(tradingtime))]
for i in range(0, len(tradingtime)):
simulated_xt_array[i] = new_fbm_simulation[int(tradingtime[i]*10)]
if GRAPHS == True:
plt.title("MMAR generated Xt")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('X(t)\n(Natural log growth since beginning)')
plt.grid(b=True)
plt.fill_between(np.arange(0, 2**K, 1) , simulated_xt_array, color="darkviolet", alpha=0.2)
plt.show()
# --- PRICES ---
if GRAPHS == True:
plt.title("MMAR generated Price")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Price level')
plt.grid(b=True)
plt.fill_between(np.arange(0, 2**K, 1) , original_price_history[0]*np.exp(simulated_xt_array), color="limegreen", alpha=0.2)
plt.show()
# --- LN CHANGES ---
if GRAPHS == True:
ln_simulated_xt_array = [0 for x in range(0, len(simulated_xt_array)-1)]
for i in range(1,len(simulated_xt_array)):
ln_simulated_xt_array[i-1] = np.log((original_price_history[0]*np.exp(simulated_xt_array[i]))/(original_price_history[0]*np.exp(simulated_xt_array[i-1])))
plt.figure(figsize=(24,5))
plt.title("Price increments")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Change\n(%)')
plt.grid(b=True)
plt.plot(ln_simulated_xt_array, color="darkviolet", linewidth=0.5)
plt.gca().set_yticklabels(['{:.0f}'.format(x*100) for x in plt.gca().get_yticks()])
plt.show()
print("The number of price increments that equal zero is: " + str(list(np.abs(ln_simulated_xt_array)).count(0)))
# --- END ---
return simulated_xt_array
In [469]:
wonderfullife = MMAR(13, H_NORWAY, lambda_NORWAY, sigma_NORWAY, df.Price_NORWAY, 0.15, True)
In [132]:
new_wonderfullife = MMAR(13, H_NORWAY, lambda_NORWAY, sigma_NORWAY, df.Price_NORWAY, 0.15, True)
In [1071]:
plt.figure(figsize=(24,5))
plt.title("Price increments")
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
axes = plt.gca()
axes.set_ylim([-0.11,0.12])
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
plt.plot(np.diff(wonderfullife), color="darkviolet", linewidth=0.5)
plt.gca().set_yticklabels(['{:.0f}'.format(x*100) for x in plt.gca().get_yticks()])
plt.show()
In [ ]:
In [65]:
# A test for creating 10,000 columns or 8,192 rows - i.e. an array with81,920,000 elements
# Apparently, everything works fine! It doesn't even take long :D
# del test_big_array
# test_big_array = [[7.0707 for x in range(0,10000)] for y in range(0,8192)]
# test_big_array = pd.DataFrame(test_big_array)
# test_big_array
In [189]:
print("We are now ready to attempt to simulate financial markets!")
In [164]:
from IPython.display import clear_output
attempt_length = 10000
attempt=pd.DataFrame([0 for x in range(0,2**13)])
start_time = datetime.now().strftime("%I:%M%p")
magnitude_NORWAY = 0.10
for i in range(0,attempt_length):
print("Having started at "+str(start_time))
if i%(attempt_length/attempt_length)==0:
print("Currently at i=" + str(i))
if i%101==0:
attempt[i] = pd.DataFrame(MMAR(13, H_NORWAY, lambda_NORWAY, sigma_NORWAY, df.Price_NORWAY, magnitude_NORWAY, False))
else:
attempt[i] = pd.DataFrame(MMAR(13, H_NORWAY, lambda_NORWAY, sigma_NORWAY, df.Price_NORWAY, magnitude_NORWAY, False))
clear_output()
print("Having started at "+str(start_time))
attempt
Out[164]:
In [173]:
# attempt.to_csv('Norway_attempt.csv')
In [154]:
plt.figure(figsize=(24,10))
plt.plot(df.Price_NORWAY, color="black", linewidth=2, label="Real Norwegian price path")
plt.legend()
for i in range(0,len(attempt.columns)):
K=13
plt.title(str(attempt_length) + " MMAR generated price paths")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Price level')
plt.grid(b=True)
plt.plot(np.arange(0, 2**K, 1) , df.Price_NORWAY[0]*np.exp(attempt[i]), color="orangered", linewidth=0.1, alpha=0.2)
plt.plot(df.Price_NORWAY, color="black", linewidth=2)
plt.show()
In [1060]:
plt.grid(b=True)
plt.title("Histogram of simulated final prices\nafter 8192 days for Norway")
plt.hist(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191]), bins=40, color="orangered")
plt.xlabel("Final price (USD/NOK)")
plt.ylabel("Frequency")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.05)), linestyle='--', color="purple", label="5th percentile")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.25)), linestyle='--', color="blue", label="25th percentile")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.5)), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.75)), linestyle='--', color="gold", label="75th percentile")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.95)), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_NORWAY.iloc[-1], linestyle='--', color="black", label="real final price", linewidth=0.7)
plt.legend()
plt.show()
In [176]:
# Here we calculate the LN differences
diff_attempt = [[0 for x in range(0,attempt_length)] for y in range (0,8192)]
for i in range(0,attempt_length):
if i(attempt_length/attempt_length) == 0:
print("i="+str(i))
for j in range (1,8192):
diff_attempt[j][i]=attempt[i][j] - attempt[i][j-1]
clear_output()
diff_attempt = pd.DataFrame(diff_attempt)
diff_attempt.to_csv('Norway_diff_attempt.csv')
In [209]:
# Here we show the original graph as a reminder
plt.figure(figsize=(24,5))
plt.title("Price increments for Norway")
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.10,0.10])
plt.plot(df.LN_NORWAY, color="red", linewidth=0.5)
plt.show()
In [210]:
# Here we show the first twenty graphs of price changes
for i in range(0,20):
plt.figure(figsize=(24,5))
plt.title("Price increments for simulation " + str(i+1))
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.10,0.10])
plt.plot(diff_attempt[i], color="darkviolet", linewidth=0.5)
plt.show()
In [179]:
print("Here we can see some statistics about the price changes in the Norwegian exchange rate.\n")
df.LN_NORWAY.describe()
Out[179]:
In [196]:
df.LN_NORWAY.std()
Out[196]:
In [1046]:
plt.grid(b=True)
plt.title("Mean price change for \n" + str(attempt_length) + " MMAR generated price paths\nfor NORWAY")
plt.hist(diff_attempt.mean(), bins=30, color='orangered')
plt.axvline(df.LN_NORWAY.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_NORWAY.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(attempt_length))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((diff_attempt.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of price changes for \n" + str(attempt_length) + " MMAR generated price paths\nfor NORWAY")
plt.hist(diff_attempt.std(), bins=30, color='orangered')
plt.axvline(df.LN_NORWAY.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_NORWAY.std(), 5)))
plt.ylabel("Frequency per " + str(attempt_length))
plt.xlabel("Standard Deviation")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our 10,000 simulations.\n")
(diff_attempt.std()).describe()
Out[1046]:
In [181]:
print("The mean is obviously very variable - we can simulate both downward and upward markets.\n\nWe can see that the standard deviation is also quite variable, but it is mostly near to what we expect.")
In [182]:
print("Here's the kurtosis for the real data for Norway.")
stats.kurtosis(df.LN_NORWAY[1:-1])
Out[182]:
In [1181]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(attempt_length) + " MMAR generated price paths\nfor Norway")
plt.hist(stats.kurtosis(diff_attempt), bins=100, color='orangered')
plt.axvline(stats.kurtosis(df.LN_NORWAY[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_NORWAY[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_diff_attempt = stats.kurtosis(diff_attempt)
plt.axvline(np.quantile(kurtosis_array_diff_attempt, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(attempt_length))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [184]:
print("But at least the kurtosis is about right!")
In [208]:
print("Below is the mean kurtosis of the price changes for every simulation.\nKeep in mind that it is overweighted by the heavy right tail.\n")
stats.kurtosis(diff_attempt).mean()
Out[208]:
In [701]:
print("The median kurtosis is about right. The real kurtosis was "+str(stats.kurtosis(df.LN_NORWAY[1:-1]))+".\n")
np.median(stats.kurtosis(diff_attempt))
Out[701]:
In [827]:
# Here we test how well the simulations fit the original data, using the Kolmogorov-Smirnov test to test the log price changes
ks_test_NORWAY = [0 for x in range(0, attempt_length)]
for i in range(0,attempt_length):
ks_test_NORWAY[i] = stats.ks_2samp(diff_attempt[i], df.LN_NORWAY).pvalue
plt.grid(b=True)
plt.hist(ks_test_NORWAY, bins=20, color='orangered')
plt.title("Histogram of Kolmogorov-Smirnov P-values\nfor " + str(attempt_length) + " MMAR simulations for Norway ")
plt.axvline(0.05, linestyle='--', color="black", label="5% significance")
plt.xlabel("P-Value")
plt.ylabel("Frequency per " + str(attempt_length))
plt.legend()
plt.show()
plt.grid(b=True)
plt.plot(sorted(ks_test_NORWAY), color='orangered')
plt.title("Sorted Kolmogorov-Smirnov P-values\nfor " + str(attempt_length) + " MMAR simulations for Norway ")
plt.axhline(0.05, linestyle='--', color="black", label="5% significance")
plt.axvline(sum(x<0.05 for x in ks_test_NORWAY), linestyle='--', color="purple", label="Total rejections = "+str(sum(x<0.05 for x in ks_test_NORWAY)))
plt.xlabel("Index ouf of " + str(attempt_length))
plt.ylabel("P-Value")
plt.legend()
plt.show()
print("Count of P<0.05 = "+ str(sum(x<0.05 for x in ks_test_NORWAY)) + " failures out of " + str(attempt_length) + ".")
In [188]:
print("We can see that most simulations fail the Kolmogorov-Smirnov test, which is a shame.\n")
print("That said, the test may be overly sensitive.\n")
print("But unfortunately, the authors don't provide any better way to test if an MMAR simulation fits the real data.")
In [ ]:
In [ ]:
In [ ]:
In [842]:
print("Now we encounter a problem: we see that the mean return for an MMAR simulation is 0.\n")
print("For stocks of rich countries, we don't expect the stock index to go down over 30 years,\nso we can't leave our price increments with an expectation of zero.")
print("\nWe must therefore introduce an upward bias into our price movement simulation.")
print("\nHere, we will multiply our simulated Xt by e to the power of the mean daily return over 30 years")
print("(in our case, the mean daily return was " + str(round(df.LN_SWEDEN.mean(), 6)*100) + "% per day.) This will bias the price to move upwards.")
print("\nBut we also don't want our prices to move in a way that's too extreme, so we will restrict our prices\nfrom ever exceeding a barrier that's 4x higher than the real price at the end of thirty years.")
print("This is admittedly an arbitray number, but it is comparable to the Nikkei 225 index, \nwhich rose to nearly 40,000 in the 80's and was still only around 10,000 over twenty years later. ")
In [839]:
def MMAR_STOCKS(K, simulated_H, simulated_lambda, simulated_sigma, original_price_history, magnitude_parameter, GRAPHS, FBMCLASS, FBM_TOP_LIMIT, FBM_BOTTOM_LIMIT, MEAN_RETURN):
# --- VARIABLES ---
# K
# adjust K depending on how many days you want to simulate (e.g. if K=13, you'll simulate 2^13=8192 days)
# simulated_H
# the Hurst parameter for the fBm process
# simulated_lambda
# the mean for the lognormal cascade
# simulated_sigma
# the variance for the lognormal cascade
# original_price_history
# the price history of the market under study (used for starting the prices from the right time!)
# magnitude_parameter
# adjust this up or down to change the range of price changes (e.g. if prices are swinging too wildly every day, then adjust this downwards)
# GRAPHS
# Boolean - either True or False - use True if you want the MMAR to simulate graphs for you
# --- MESSAGE ---
if GRAPHS == True:
print("Performing an MMAR simulation with parameters:\n\nH = " + str(simulated_H) + "\nlambda = " + str(simulated_lambda) + "\nsigma = " + str(simulated_sigma) + "\nfBm magnitude = " + str(magnitude_parameter)+ "\n")
# --- CASCADE ---
new_cascade = list(np.array(lognormal_cascade(k=K, v=1, ln_lambda = simulated_lambda, ln_theta = simulated_sigma)).flat)
if GRAPHS == True:
plt.figure(figsize=(24,2))
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.title("Binomial Cascade")
plt.xlabel("Conventional time\n(days)")
plt.ylabel('"Mass"')
plt.plot(new_cascade, color="crimson", linewidth=0.5)
plt.show()
# --- TRADING TIME ---
tradingtime = 2**K*np.cumsum(new_cascade)/sum(new_cascade)
if GRAPHS == True:
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.yticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.title("Trading time\n(normalized)")
plt.xlabel("Conventional time\n(days)")
plt.ylabel('"Trading time"\n(\u03B8 (t), normalized)')
plt.plot(tradingtime, color="orangered")
# --- FBM (Fractional Brownian Motion) ---
# Here we are adding the complicating factor. It is a recursive method to keep simulating fBm's until we get one that gives us a realistic price movement.
def fbm_attempt(fbmclass, fbm_top_limit, fbm_bottom_limit, mean_return, original_prices):
# --- FBM (Fractional Brownian Motion) ---
simulated_fbm = fbmclass.fbm()
xt_top_limit = np.log(fbm_top_limit/original_prices[0])
xt_bottom_limit = np.log(fbm_bottom_limit/original_prices[0])
# --- MMAR XT's (temporary)---
TEMPORARY_simulated_xt_array = [simulated_fbm[int(tradingtime[x]*10)]+x*mean_return for x in range(0, len(tradingtime))] # Here, since it's stocks, we also add a multiple of the mean return for every Xt value
if all(XT < xt_top_limit for XT in TEMPORARY_simulated_xt_array) and all(XT > xt_bottom_limit for XT in TEMPORARY_simulated_xt_array):
return TEMPORARY_simulated_xt_array
# --- Graphs for FBM ---
if GRAPHS == True:
plt.figure(figsize=(24,2))
plt.xticks(np.arange(0, 10*2**(K)+1, 10*2**(K-3)))
plt.title("Fractional Brownian Motion")
plt.xlabel("t")
plt.ylabel('fBm (t)')
plt.plot(simulated_fbm, color="orange")
plt.show()
else:
return fbm_attempt(fbmclass, fbm_top_limit, fbm_bottom_limit, mean_return, original_prices)
# --- MMAR XT's ---
simulated_xt_array = fbm_attempt(FBMCLASS, FBM_TOP_LIMIT, FBM_BOTTOM_LIMIT, MEAN_RETURN, original_price_history)
# --- Graphs for MMAR XT's ---
if GRAPHS == True:
plt.title("MMAR generated Xt")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('X(t)\n(Natural log growth since beginning)')
plt.grid(b=True)
plt.fill_between(np.arange(0, 2**K, 1) , simulated_xt_array, color="darkviolet", alpha=0.2)
plt.show()
# --- PRICES ---
if GRAPHS == True:
plt.title("MMAR generated Price")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Price level')
plt.grid(b=True)
plt.fill_between(np.arange(0, 2**K, 1) , original_price_history[0]*np.exp(simulated_xt_array), color="limegreen", alpha=0.2)
plt.show()
# --- LN CHANGES ---
if GRAPHS == True:
ln_simulated_xt_array = [0 for x in range(0, len(simulated_xt_array)-1)]
for i in range(1,len(simulated_xt_array)):
ln_simulated_xt_array[i-1] = np.log((original_price_history[0]*np.exp(simulated_xt_array[i]))/(original_price_history[0]*np.exp(simulated_xt_array[i-1])))
plt.figure(figsize=(24,5))
plt.title("Price increments")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Change\n(%)')
plt.grid(b=True)
plt.plot(ln_simulated_xt_array, color="darkviolet", linewidth=0.5)
plt.gca().set_yticklabels(['{:.0f}'.format(x*100) for x in plt.gca().get_yticks()])
plt.show()
print("The number of price increments that equal zero is: " + str(list(np.abs(ln_simulated_xt_array)).count(0)))
# --- END ---
return simulated_xt_array
In [992]:
from IPython.display import clear_output
attempt2_length = 10000
attempt2=pd.DataFrame([0 for x in range(0,2**13)])
magnitude_SWEDEN = 3.24
FBM_class_SWEDEN = FBM(n = 10*2**13+1, hurst = H_SWEDEN, length = magnitude_SWEDEN, method='daviesharte')
start_time = datetime.now().strftime("%I:%M%p")
for i in range(0,attempt2_length):
print("Having started at "+str(start_time))
if i%(attempt2_length/attempt2_length)==0:
print("Currently at i=" + str(i) + " out of " + str(attempt2_length))
if i%101==0:
attempt2[i] = pd.DataFrame(MMAR_STOCKS(13, H_SWEDEN, lambda_SWEDEN, sigma_SWEDEN, df.Price_SWEDEN, magnitude_SWEDEN, GRAPHS=False, FBMCLASS=FBM_class_SWEDEN, FBM_TOP_LIMIT=5000, FBM_BOTTOM_LIMIT=100, MEAN_RETURN=df.LN_SWEDEN.mean()))
else:
attempt2[i] = pd.DataFrame(MMAR_STOCKS(13, H_SWEDEN, lambda_SWEDEN, sigma_SWEDEN, df.Price_SWEDEN, magnitude_SWEDEN, GRAPHS=False, FBMCLASS=FBM_class_SWEDEN, FBM_TOP_LIMIT=5000, FBM_BOTTOM_LIMIT=100, MEAN_RETURN=df.LN_SWEDEN.mean()))
clear_output()
print("Having started at "+str(start_time))
attempt2
Out[992]:
In [993]:
attempt2.to_csv('Sweden_attempt.csv')
In [155]:
plt.figure(figsize=(24,10))
plt.plot(df.Price_SWEDEN, color="black", linewidth=2, label="Real Swedish price path")
plt.legend()
for i in range(0,len(attempt2.columns)):
K=13
plt.title(str(attempt2_length) + " MMAR generated price paths")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Price level')
plt.grid(b=True)
plt.plot(np.arange(0, 2**K, 1) , df.Price_SWEDEN[0]*np.exp(attempt2[i]), color="goldenrod", linewidth=0.1, alpha=0.2)
plt.plot(df.Price_SWEDEN, color="black", linewidth=2)
plt.show()
In [1052]:
plt.grid(b=True)
plt.title("Histogram of simulated final prices\nafter 8192 days for Sweden")
plt.hist(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191]), bins=40, color="goldenrod")
plt.xlabel("Final price (SEK)")
plt.ylabel("Frequency")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.05)), linestyle='--', color="purple", label="5th percentile")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.25)), linestyle='--', color="blue", label="25th percentile")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.5)), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.75)), linestyle='--', color="gold", label="75th percentile")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.95)), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_SWEDEN.iloc[-1], linestyle='--', color="black", label="real final price")
plt.legend()
plt.show()
In [996]:
# Here we calculate the LN differences
diff_attempt2 = [[0 for x in range(0,attempt2_length)] for y in range (0,8192)]
for i in range(0,attempt2_length):
if i%(attempt2_length/attempt2_length) == 0:
print("i="+str(i))
for j in range (1,8192):
diff_attempt2[j][i]=attempt2[i][j] - attempt2[i][j-1]
clear_output()
diff_attempt2= pd.DataFrame(diff_attempt2)
diff_attempt2.to_csv('Sweden_diff_attempt.csv')
In [997]:
# Here we show the original graph as a reminder
plt.figure(figsize=(24,5))
plt.title("Price increments for Sweden")
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.175,0.175])
plt.plot(df.LN_SWEDEN, color="gold", linewidth=0.5)
plt.show()
In [998]:
# Here we show the first twenty graphs of price changes
for i in range(0,20):
plt.figure(figsize=(24,5))
plt.title("Price increments for simulation " + str(i+1))
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.175,0.175])
plt.plot(diff_attempt2[i], color="orange", linewidth=0.5)
plt.show()
In [999]:
print("Here we can see some statistics about the price changes in the Swedish stock market.\n")
df.LN_SWEDEN.describe()
Out[999]:
In [1000]:
df.LN_SWEDEN.std()
Out[1000]:
In [1036]:
plt.grid(b=True)
plt.title("Mean price change for \n" + str(attempt2_length) + " MMAR generated price paths\nfor SWEDEN")
plt.hist(diff_attempt2.mean(), bins=30, color='goldenrod')
plt.axvline(df.LN_SWEDEN.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_SWEDEN.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(attempt2_length))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((diff_attempt2.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of price changes for \n" + str(attempt2_length) + " MMAR generated price paths\nfor SWEDEN")
plt.hist(diff_attempt2.std(), bins=30, color='goldenrod')
plt.axvline(df.LN_SWEDEN.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_SWEDEN.std(), 4)))
plt.ylabel("Frequency per " + str(attempt2_length))
plt.xlabel("Standard Deviation")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our " + str(attempt2_length) + " simulations.\n")
(diff_attempt2.std()).describe()
Out[1036]:
In [1002]:
print("The mean is obviously very variable - we can simulate both downward and upward markets.\n\nWe can see that the standard deviation is also quite variable, but it is mostly near to what we expect.")
In [1003]:
print("Here's the kurtosis for the real data for Sweden.")
stats.kurtosis(df.LN_SWEDEN[1:-1])
Out[1003]:
In [1180]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(attempt2_length) + " MMAR generated price paths\nfor Sweden")
plt.hist(stats.kurtosis(diff_attempt2), bins=100, color='goldenrod')
plt.axvline(stats.kurtosis(df.LN_SWEDEN[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_SWEDEN[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_diff_attempt2 = stats.kurtosis(diff_attempt2)
plt.axvline(np.quantile(kurtosis_array_diff_attempt2, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt2,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt2,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt2,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt2,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(attempt2_length))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [1005]:
print("But at least the kurtosis is about right!")
In [1006]:
print("Below is the mean kurtosis of the price changes for every simulation.\nKeep in mind that it is overweighted by the heavy right tail.\n")
stats.kurtosis(diff_attempt2).mean()
Out[1006]:
In [1007]:
print("The median kurtosis is a little less than the real, which was "+str(stats.kurtosis(df.LN_SWEDEN[1:-1]))+".\n")
np.median(stats.kurtosis(diff_attempt2))
Out[1007]:
In [1008]:
# Here we test how well the simulations fit the original data, using the Kolmogorov-Smirnov test to test the log price changes
ks_test_SWEDEN = [0 for x in range(0, attempt2_length)]
for i in range(0,attempt2_length):
ks_test_SWEDEN[i] = stats.ks_2samp(diff_attempt2[i], df.LN_SWEDEN).pvalue
plt.grid(b=True)
plt.hist(ks_test_SWEDEN, bins=20, color='goldenrod')
plt.title("Histogram of Kolmogorov-Smirnov P-values\nfor " + str(attempt2_length) + " MMAR simulations for Sweden ")
plt.axvline(0.05, linestyle='--', color="black", label="5% significance")
plt.xlabel("P-Value")
plt.ylabel("Frequency per " + str(attempt2_length))
plt.legend()
plt.show()
plt.grid(b=True)
plt.plot(sorted(ks_test_SWEDEN), color='goldenrod')
plt.title("Sorted Kolmogorov-Smirnov P-values\nfor " + str(attempt2_length) + " MMAR simulations for Sweden ")
plt.axhline(0.05, linestyle='--', color="black", label="5% significance")
plt.axvline(sum(x<0.05 for x in ks_test_SWEDEN), linestyle='--', color="orange", label="Total rejections = "+str(sum(x<0.05 for x in ks_test_SWEDEN)))
plt.xlabel("Index ouf of " + str(attempt2_length))
plt.ylabel("P-Value")
plt.legend()
plt.show()
print("Count of P<0.05 = "+ str(sum(x<0.05 for x in ks_test_SWEDEN)) + " failures out of " + str(attempt2_length) + ".")
In [1009]:
print("We can see that most simulations fail the Kolmogorov-Smirnov test, which is a shame.\n")
print("That said, the test may be overly sensitive.\n")
print("But unfortunately, the authors don't provide any better way to test if an MMAR simulation fits the real data.")
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [1010]:
print("In the same way that we expect stocks to go downwards, we may expect lending rates to move in a certain direction.")
print("\nHere, we reuse our modified MMAR for stocks to simulate bond yields going down over thirty years.")
print("In the same way as for stocks, we can modify our Xt's by adding the mean change from the real yield changes over thirty years.")
In [1011]:
from IPython.display import clear_output
attempt3_length = 10000
attempt3=pd.DataFrame([0 for x in range(0,2**13)])
magnitude_LIBOR = 4
FBM_class_LIBOR = FBM(n = 10*2**13+1, hurst = H_LIBOR, length = magnitude_LIBOR, method='daviesharte')
start_time = datetime.now().strftime("%I:%M%p")
for i in range(0,attempt3_length):
print("Having started at "+str(start_time))
if i%(attempt3_length/attempt3_length)==0:
print("Currently at i=" + str(i) + " out of " + str(attempt3_length))
if i%101==0:
attempt3[i] = pd.DataFrame(MMAR_STOCKS(13, H_LIBOR, lambda_LIBOR, sigma_LIBOR, df.Price_LIBOR, magnitude_LIBOR, GRAPHS=False, FBMCLASS=FBM_class_LIBOR, FBM_TOP_LIMIT=20, FBM_BOTTOM_LIMIT=0.1, MEAN_RETURN=df.LN_LIBOR.mean()))
else:
attempt3[i] = pd.DataFrame(MMAR_STOCKS(13, H_LIBOR, lambda_LIBOR, sigma_LIBOR, df.Price_LIBOR, magnitude_LIBOR, GRAPHS=False, FBMCLASS=FBM_class_LIBOR, FBM_TOP_LIMIT=20, FBM_BOTTOM_LIMIT=0.1, MEAN_RETURN=df.LN_LIBOR.mean()))
clear_output()
print("Having started at "+str(start_time))
attempt3
Out[1011]:
In [1012]:
attempt3.to_csv('LIBOR_attempt.csv')
In [149]:
plt.figure(figsize=(24,10))
plt.plot(df.Price_LIBOR, color="black", linewidth=2, label="Real LIBOR price path")
plt.legend()
for i in range(0,len(attempt3.columns)):
K=13
plt.title(str(attempt3_length) + " MMAR generated price paths")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Price level')
plt.grid(b=True)
plt.plot(np.arange(0, 2**K, 1) , df.Price_LIBOR[0]*np.exp(attempt3[i]), color="cornflowerblue", linewidth=0.1, alpha=0.2)
plt.plot(df.Price_LIBOR, color="black", linewidth=2)
plt.show()
In [1055]:
plt.grid(b=True)
plt.title("Histogram of simulated final yields\nafter 8192 days for LIBOR")
plt.hist(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191]), bins=40, color="cornflowerblue")
plt.xlabel("Final yield (%)")
plt.ylabel("Frequency")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.05)), linestyle='--', color="purple", label="5th percentile")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.25)), linestyle='--', color="blue", label="25th percentile")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.5)), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.75)), linestyle='--', color="gold", label="75th percentile")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.95)), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_LIBOR.iloc[-1], linestyle='--', color="black", label="real final price")
plt.legend()
plt.show()
In [1014]:
# Here we calculate the LN differences
diff_attempt3 = [[0 for x in range(0,attempt3_length)] for y in range (0,8192)]
for i in range(0,attempt3_length):
if i%(attempt3_length/attempt3_length) == 0:
print("i="+str(i))
for j in range (1,8192):
diff_attempt3[j][i]=attempt3[i][j] - attempt3[i][j-1]
clear_output()
diff_attempt3= pd.DataFrame(diff_attempt3)
diff_attempt3.to_csv('LIBOR_diff_attempt.csv')
In [1015]:
# Here we show the original graph as a reminder
plt.figure(figsize=(24,5))
plt.title("Price increments for LIBOR")
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.25,0.25])
plt.plot(df.LN_LIBOR, color="blue", linewidth=0.5)
plt.show()
In [1016]:
# Here we show the first twenty graphs of price changes
for i in range(0,20):
plt.figure(figsize=(24,5))
plt.title("Price increments for simulation " + str(i+1))
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.25,0.25])
plt.plot(diff_attempt3[i], color="forestgreen", linewidth=0.5)
plt.show()
In [1017]:
print("Here we can see some statistics about the price changes in the LIBOR rate in British pounds.\n")
df.LN_LIBOR.describe()
Out[1017]:
In [1018]:
df.LN_LIBOR.std()
Out[1018]:
In [1043]:
plt.grid(b=True)
plt.title("Mean yield change for \n" + str(attempt3_length) + " MMAR generated price paths\nfor 12-month LIBOR")
plt.hist(diff_attempt3.mean(), bins=30, color='cornflowerblue')
plt.axvline(df.LN_LIBOR.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_LIBOR.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(attempt3_length))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((diff_attempt3.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of yield changes for \n" + str(attempt3_length) + " MMAR generated price paths\nfor 12-month LIBOR")
plt.hist(diff_attempt3.std(), bins=30, color='cornflowerblue')
plt.axvline(df.LN_LIBOR.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_LIBOR.std(), 5)))
plt.ylabel("Frequency per " + str(attempt3_length))
plt.xlabel("Standard Deviation")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our " + str(attempt3_length) + " simulations.\n")
print((diff_attempt3.std()).describe())
In [1020]:
print("The mean is obviously very variable - we can simulate both downward and upward markets.\n\nWe can see that the standard deviation is also quite variable, but it is mostly near to what we expect.")
In [1021]:
print("Here's the kurtosis for the real data for LIBOR.")
stats.kurtosis(df.LN_LIBOR[1:-1])
Out[1021]:
In [1179]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of yield changes for \n" + str(attempt3_length) + " MMAR generated price paths\nfor 12-month LIBOR")
plt.hist(stats.kurtosis(diff_attempt3), bins=100, color='cornflowerblue')
plt.axvline(stats.kurtosis(df.LN_LIBOR[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_LIBOR[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_diff_attempt3 = stats.kurtosis(diff_attempt3)
plt.axvline(np.quantile(kurtosis_array_diff_attempt3, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt3,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt3,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt3,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_diff_attempt3,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(attempt3_length))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [1023]:
print("For LIBOR, the kurtosis is kind of all over the place!")
In [1024]:
print("Below is the mean kurtosis of the price changes for every simulation.\nKeep in mind that it is overweighted by the heavy right tail.\n")
stats.kurtosis(diff_attempt3).mean()
Out[1024]:
In [1025]:
print("The median kurtosis is quite far from the real value, which was "+str(stats.kurtosis(df.LN_LIBOR[1:-1]))+".\n")
np.median(stats.kurtosis(diff_attempt3))
Out[1025]:
In [1067]:
# Here we test how well the simulations fit the original data, using the Kolmogorov-Smirnov test to test the log price changes
ks_test_LIBOR = [0 for x in range(0, attempt3_length)]
for i in range(0,attempt3_length):
ks_test_LIBOR[i] = stats.ks_2samp(diff_attempt3[i], df.LN_LIBOR).pvalue
plt.grid(b=True)
plt.hist(ks_test_LIBOR, bins=20, color='cornflowerblue')
plt.title("Histogram of Kolmogorov-Smirnov P-values\nfor " + str(attempt3_length) + " MMAR simulations for LIBOR ")
plt.axvline(0.05, linestyle='--', color="black", label="5% significance")
plt.xlabel("P-Value")
plt.ylabel("Frequency per " + str(attempt3_length))
plt.legend()
plt.show()
plt.grid(b=True)
plt.plot(sorted(ks_test_LIBOR), color='cornflowerblue')
plt.title("Sorted Kolmogorov-Smirnov P-values\nfor " + str(attempt3_length) + " MMAR simulations for LIBOR ")
plt.axhline(0.05, linestyle='--', color="black", label="5% significance")
plt.axvline(sum(x<0.05 for x in ks_test_LIBOR), linestyle='--', color="forestgreen", label="Total rejections = "+str(sum(x<0.05 for x in ks_test_LIBOR)))
plt.xlabel("Index ouf of " + str(attempt3_length))
plt.ylabel("P-Value")
plt.legend()
plt.show()
print("Count of P<0.05 = "+ str(sum(x<0.05 for x in ks_test_LIBOR)) + " failures out of " + str(attempt3_length) + ".")
In [1027]:
print("We can see that most simulations fail the Kolmogorov-Smirnov test, which is a shame.\n")
print("That said, the test may be overly sensitive.\n")
print("But unfortunately, the authors don't provide any better way to test if an MMAR simulation fits the real data.")
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [965]:
K=13
for i in range(0,1000):
if i%500 == 0:
print("currently at i=" + str(i))
cherie = np.cumsum(lognormal_cascade(k=K, v=1, ln_lambda=lambda_NORWAY, ln_theta=sigma_NORWAY))
cherie = 2**(K)*cherie/cherie[-1]
plt.plot(cherie, linewidth=0.1, color="orangered")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.yticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.title("1,000 simulations of trading time\nfor Norway")
plt.xlabel("Conventional time\n(days)")
plt.ylabel('"Trading time"\n(\u03B8 (t), normalized)')
plt.grid(b=True)
plt.show()
for i in range(0,1000):
if i%500 == 0:
print("currently at i=" + str(i))
sweetie = np.cumsum(lognormal_cascade(k=K, v=1, ln_lambda=lambda_SWEDEN, ln_theta=sigma_SWEDEN))
sweetie = 2**(K)*sweetie/sweetie[-1]
plt.plot(sweetie, linewidth=0.1, color="goldenrod")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.yticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.title("1,000 simulations of trading time\nfor Sweden")
plt.xlabel("Conventional time\n(days)")
plt.ylabel('"Trading time"\n(\u03B8 (t), normalized)')
plt.grid(b=True)
plt.show()
for i in range(0,1000):
if i%500 == 0:
print("currently at i=" + str(i))
libbie = np.cumsum(lognormal_cascade(k=K, v=1, ln_lambda=lambda_LIBOR, ln_theta=sigma_LIBOR))
libbie = 2**(K)*libbie/libbie[-1]
plt.plot(libbie, linewidth=0.1, color="cornflowerblue")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.yticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.title("1,000 simulations of trading time\nfor LIBOR")
plt.xlabel("Conventional time\n(days)")
plt.ylabel('"Trading time"\n(\u03B8 (t), normalized)')
plt.grid(b=True)
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [1161]:
Gaussian_length = 10000
Gaussian_Norway = pd.DataFrame()
start_time = datetime.now().strftime("%I:%M%p")
for i in range(0,Gaussian_length):
print("Having started at "+str(start_time))
print("i="+str(i))
Gaussian_Norway[i] = [0.0 for x in range(0,2**13)]
Gaussian_Norway[i][0] = df.Price_NORWAY[0]
for j in range(1,len(Gaussian_Norway)):
Gaussian_Norway[i][j] = Gaussian_Norway[i][j-1]*np.exp(np.random.normal(0,df.LN_NORWAY.std()))
clear_output()
print("Having started at "+str(start_time))
Gaussian_Norway.to_csv('Gaussian_Norway.csv')
In [1162]:
plt.figure(figsize=(24,10))
for i in range(0,Gaussian_length):
K=13
plt.title(str(Gaussian_length) + " Gaussian generated price paths")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Price level')
plt.grid(b=True)
plt.plot(Gaussian_Norway[i], color="orangered", linewidth=0.1, alpha=0.2)
plt.show()
In [1163]:
# Here we calculate the LN differences
Gaussian_diff_attempt = [[0 for x in range(0,Gaussian_length)] for y in range (0,8192)]
for i in range(0,Gaussian_length):
if i%(Gaussian_length/Gaussian_length) == 0:
print("i="+str(i))
Gaussian_diff_attempt[i]=np.log(Gaussian_Norway[i]).diff()
clear_output()
Gaussian_diff_attempt= pd.DataFrame(Gaussian_diff_attempt)
Gaussian_diff_attempt.to_csv('Gaussian_diff_Norway.csv')
In [76]:
for i in range(0,5):
plt.figure(figsize=(24,5))
plt.title("Price increments for simulation " + str(i+1))
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.10,0.10])
plt.plot(Gaussian_diff_attempt[i], color="darkviolet", linewidth=0.5)
plt.show()
In [1174]:
plt.grid(b=True)
plt.title("Histogram of Gaussian simulated final prices\nafter 8192 days for Norway")
plt.hist(Gaussian_Norway.iloc[-1], bins=40, color="orangered")
plt.xlabel("Final price (USD/NOK)")
plt.ylabel("Frequency")
plt.axvline(Gaussian_Norway.iloc[-1].quantile(0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(Gaussian_Norway.iloc[-1].quantile(0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(Gaussian_Norway.iloc[-1].quantile(0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(Gaussian_Norway.iloc[-1].quantile(0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(Gaussian_Norway.iloc[-1].quantile(0.95), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_NORWAY.iloc[-1], linestyle='--', color="black", label="real final price", linewidth=0.7)
plt.legend()
plt.show()
In [179]:
print("Here we can see some statistics about the price changes in the Norwegian exchange rate.\n")
df.LN_NORWAY.describe()
Out[179]:
In [196]:
df.LN_NORWAY.std()
Out[196]:
In [67]:
attempt = pd.read_csv("~sib-ros/Desktop/Simulation collection/Norway_attempt.csv", index_col=0)
attempt.columns = attempt.columns.astype(int)
diff_attempt = pd.read_csv("~sib-ros/Desktop/Simulation collection/Norway_diff_attempt.csv", index_col=0)
diff_attempt.columns = diff_attempt.columns.astype(int)
Gaussian_Norway = pd.read_csv("~sib-ros/Desktop/Simulation collection/Gaussian_Norway.csv", index_col=0)
Gaussian_Norway.columns = attempt.columns.astype(int)
Gaussian_diff_attempt = pd.read_csv("~sib-ros/Desktop/Simulation collection/Gaussian_diff_Norway.csv", index_col=0)
Gaussian_diff_attempt.columns = diff_attempt.columns.astype(int)
In [185]:
plt.grid(b=True)
plt.title("Mean price change for \n" + str(Gaussian_length) + " Gaussian generated price paths\nfor NORWAY")
plt.hist(Gaussian_diff_attempt.mean(), bins=30, color='orangered')
plt.axvline(df.LN_NORWAY.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_NORWAY.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(Gaussian_length))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((Gaussian_diff_attempt.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of price changes for \n" + str(Gaussian_length) + " Gaussian generated price paths\nfor NORWAY")
plt.hist(Gaussian_diff_attempt.std(), bins=30, color='orangered')
plt.axvline(df.LN_NORWAY.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_NORWAY.std(), 5)))
plt.ylabel("Frequency per " + str(Gaussian_length))
plt.xlabel("Standard Deviation")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our 10,000 simulations.\n")
(Gaussian_diff_attempt.std()).describe()
Out[185]:
In [181]:
print("The mean is obviously very variable - we can simulate both downward and upward markets.\n\nWe can see that the standard deviation is also quite variable, but it is mostly near to what we expect.")
In [182]:
print("Here's the kurtosis for the real data for Norway.")
stats.kurtosis(df.LN_NORWAY[1:-1])
Out[182]:
In [88]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length) + " Gaussian generated price paths\nfor Norway")
plt.hist(stats.kurtosis(Gaussian_diff_attempt), bins=100, color='orangered')
plt.axvline(stats.kurtosis(df.LN_NORWAY[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_NORWAY[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt = stats.kurtosis(Gaussian_diff_attempt)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [87]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length) + " Gaussian generated price paths\nfor Norway")
plt.hist(stats.kurtosis(Gaussian_diff_attempt), bins=100, color='orangered')
# plt.axvline(stats.kurtosis(df.LN_NORWAY[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_NORWAY[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt = stats.kurtosis(Gaussian_diff_attempt)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [100]:
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length) + " Gaussian generated price paths\nfor Norway")
plt.hist(stats.kurtosis(Gaussian_diff_attempt), bins=100, color='orangered')
# plt.axvline(stats.kurtosis(df.LN_NORWAY[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_NORWAY[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt = stats.kurtosis(Gaussian_diff_attempt)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [184]:
print("But at least the kurtosis is about right!")
In [208]:
print("Below is the mean kurtosis of the price changes for every simulation.\nKeep in mind that it is overweighted by the heavy right tail.\n")
stats.kurtosis(diff_attempt).mean()
Out[208]:
In [701]:
print("The median kurtosis is about right. The real kurtosis was "+str(stats.kurtosis(df.LN_NORWAY[1:-1]))+".\n")
np.median(stats.kurtosis(diff_attempt))
Out[701]:
In [1182]:
# Here we test how well the simulations fit the original data, using the Kolmogorov-Smirnov test to test the log price changes
ks_test_Gaussian_NORWAY = [0 for x in range(0, Gaussian_length)]
for i in range(0,Gaussian_length):
ks_test_Gaussian_NORWAY[i] = stats.ks_2samp(Gaussian_diff_attempt[i], df.LN_NORWAY).pvalue
plt.grid(b=True)
plt.hist(ks_test_Gaussian_NORWAY, bins=20, color='orangered')
plt.title("Histogram of Kolmogorov-Smirnov P-values\nfor " + str(Gaussian_length) + " Gaussian simulations for Norway ")
plt.axvline(0.05, linestyle='--', color="black", label="5% significance")
plt.xlabel("P-Value")
plt.ylabel("Frequency per " + str(Gaussian_length))
plt.legend()
plt.show()
plt.grid(b=True)
plt.plot(sorted(ks_test_Gaussian_NORWAY), color='orangered')
plt.title("Sorted Kolmogorov-Smirnov P-values\nfor " + str(Gaussian_length) + " Gaussian simulations for Norway ")
plt.axhline(0.05, linestyle='--', color="black", label="5% significance")
plt.axvline(sum(x<0.05 for x in ks_test_Gaussian_NORWAY), linestyle='--', color="purple", label="Total rejections = "+str(sum(x<0.05 for x in ks_test_Gaussian_NORWAY)))
plt.xlabel("Index ouf of " + str(Gaussian_length))
plt.ylabel("P-Value")
plt.legend()
plt.show()
print("Count of P<0.05 = "+ str(sum(x<0.05 for x in ks_test_Gaussian_NORWAY)) + " failures out of " + str(Gaussian_length) + ".")
In [1183]:
print("We can see that all Gaussian simulations fail the Kolmogorov-Smirnov test.\n")
In [1216]:
print("Like last time, we want to add some restrictions to our simulations, to avoid 'absurd' results.")
print("\nWe use the same restrictions as last time, namely:\n")
print("For Sweden: 100 < P(t) < 5000 for all t.")
print("For LIBOR: 0.1 < P(t) < 20 for all t.")
In [1197]:
def gaussian_attempt(top_limit, bottom_limit, mean_return, std_return, original_price, number_of_simulated_days):
TEMPORARY_array = [0.0 for x in range(number_of_simulated_days)]
for j in range(0,number_of_simulated_days):
if j == 0:
TEMPORARY_array[j] = original_price
else:
Z = np.random.normal(mean_return, std_return)
TEMPORARY_array[j] = TEMPORARY_array[j-1]*np.exp(Z)
if (all(Pt < top_limit for Pt in TEMPORARY_array) and all(Pt > bottom_limit for Pt in TEMPORARY_array))==True:
return TEMPORARY_array
else:
return gaussian_attempt(top_limit, bottom_limit, mean_return, std_return, original_price, number_of_simulated_days)
In [1241]:
print("We can now beging generating simple, 'random walk'-style market simulations.")
In [1242]:
Gaussian_length2 = 10000
Gaussian_Sweden = pd.DataFrame()
Sweden_top_limit = 5000
Sweden_bottom_limit = 100
number_of_simulated_days = 2**13+1
start_time = datetime.now().strftime("%I:%M%p")
for i in range(0,Gaussian_length2):
print("Having started at "+str(start_time))
print("i="+str(i))
Gaussian_Sweden[i] = gaussian_attempt(Sweden_top_limit, Sweden_bottom_limit, df.LN_SWEDEN.mean(), df.LN_SWEDEN.std(), df.Price_SWEDEN[0], number_of_simulated_days)
clear_output()
print("Having started at "+str(start_time))
Gaussian_Sweden.to_csv('Gaussian_Sweden.csv')
In [1243]:
plt.figure(figsize=(24,10))
for i in range(0,Gaussian_length2):
K=13
plt.title(str(Gaussian_length2) + " Gaussian generated price paths")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Price level')
plt.grid(b=True)
plt.plot(Gaussian_Sweden[i], color="goldenrod", linewidth=0.1, alpha=0.2)
plt.show()
In [1244]:
# Here we calculate the LN differences
Gaussian_diff_attempt2 = [np.log(Gaussian_Sweden[x]).diff() for x in range(0,Gaussian_length2)]
Gaussian_diff_attempt2 = pd.DataFrame(Gaussian_diff_attempt2).transpose()
Gaussian_diff_attempt2 = Gaussian_diff_attempt2[:-1]
Gaussian_diff_attempt2
Gaussian_diff_attempt2.to_csv('Gaussian_diff_Sweden.csv')
In [70]:
for i in range(0,5):
plt.figure(figsize=(24,5))
plt.title("Price increments for simulation " + str(i+1))
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.175,0.175])
plt.plot(Gaussian_diff_attempt2[i], color="orange", linewidth=0.5)
plt.show()
In [93]:
plt.grid(b=True)
plt.title("Histogram of Gaussian simulated final prices\nafter 8192 days for Sweden")
plt.hist(Gaussian_Sweden.iloc[-1], bins=40, color="goldenrod")
plt.xlabel("Final price (SEK)")
plt.ylabel("Frequency")
plt.axvline(Gaussian_Sweden.iloc[-1].quantile(0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(Gaussian_Sweden.iloc[-1].quantile(0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(Gaussian_Sweden.iloc[-1].quantile(0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(Gaussian_Sweden.iloc[-1].quantile(0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(Gaussian_Sweden.iloc[-1].quantile(0.95), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_SWEDEN.iloc[-1], linestyle='--', color="black", label="real final price", linewidth=0.7)
plt.legend()
plt.show()
In [77]:
plt.grid(b=True)
plt.title("Mean price change for \n" + str(Gaussian_length2) + " Gaussian generated price paths\nfor SWEDEN")
plt.hist(Gaussian_diff_attempt2.mean(), bins=30, color='goldenrod')
plt.axvline(df.LN_SWEDEN.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_SWEDEN.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(Gaussian_length2))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((Gaussian_diff_attempt2.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of price changes for \n" + str(Gaussian_length2) + " Gaussian generated price paths\nfor SWEDEN")
plt.hist(Gaussian_diff_attempt2.std(), bins=30, color='goldenrod')
plt.axvline(df.LN_SWEDEN.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_SWEDEN.std(), 5)))
plt.ylabel("Frequency per " + str(Gaussian_length2))
plt.xlabel("Standard Deviation")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our 10,000 simulations.\n")
(Gaussian_diff_attempt2.std()).describe()
Out[77]:
In [181]:
print("The mean is obviously very variable - we can simulate both downward and upward markets.\n\nWe can see that the standard deviation is also quite variable, but it is mostly near to what we expect.")
In [182]:
print("Here's the kurtosis for the real data for Norway.")
stats.kurtosis(df.LN_NORWAY[1:-1])
Out[182]:
In [79]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length2) + " Gaussian generated price paths\nfor Sweden")
plt.hist(stats.kurtosis(Gaussian_diff_attempt2), bins=100, color='goldenrod')
plt.axvline(stats.kurtosis(df.LN_SWEDEN[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_SWEDEN[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt2 = stats.kurtosis(Gaussian_diff_attempt2)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length2))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [91]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length2) + " Gaussian generated price paths\nfor Sweden")
plt.hist(stats.kurtosis(Gaussian_diff_attempt2), bins=100, color='goldenrod')
# plt.axvline(stats.kurtosis(df.LN_SWEDEN[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_SWEDEN[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt2 = stats.kurtosis(Gaussian_diff_attempt2)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length2))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [101]:
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length2) + " Gaussian generated price paths\nfor Sweden")
plt.hist(stats.kurtosis(Gaussian_diff_attempt2), bins=100, color='goldenrod')
# plt.axvline(stats.kurtosis(df.LN_SWEDEN[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_SWEDEN[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt2 = stats.kurtosis(Gaussian_diff_attempt2)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length2))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [80]:
# Here we test how well the simulations fit the original data, using the Kolmogorov-Smirnov test to test the log price changes
ks_test_Gaussian_SWEDEN = [0 for x in range(0, Gaussian_length2)]
for i in range(0,Gaussian_length2):
ks_test_Gaussian_SWEDEN[i] = stats.ks_2samp(Gaussian_diff_attempt2[i], df.LN_SWEDEN).pvalue
plt.grid(b=True)
plt.hist(ks_test_Gaussian_SWEDEN, bins=20, color='goldenrod')
plt.title("Histogram of Kolmogorov-Smirnov P-values\nfor " + str(Gaussian_length2) + " Gaussian simulations for Sweden")
plt.axvline(0.05, linestyle='--', color="black", label="5% significance")
plt.xlabel("P-Value")
plt.ylabel("Frequency per " + str(Gaussian_length2))
plt.legend()
plt.show()
plt.grid(b=True)
plt.plot(sorted(ks_test_Gaussian_SWEDEN), color='goldenrod')
plt.title("Sorted Kolmogorov-Smirnov P-values\nfor " + str(Gaussian_length2) + " Gaussian simulations for Sweden")
plt.axhline(0.05, linestyle='--', color="black", label="5% significance")
plt.axvline(sum(x<0.05 for x in ks_test_Gaussian_SWEDEN), linestyle='--', color="orange", label="Total rejections = "+str(sum(x<0.05 for x in ks_test_Gaussian_SWEDEN)))
plt.xlabel("Index ouf of " + str(Gaussian_length2))
plt.ylabel("P-Value")
plt.legend()
plt.show()
print("Count of P<0.05 = "+ str(sum(x<0.05 for x in ks_test_Gaussian_SWEDEN)) + " failures out of " + str(Gaussian_length2) + ".")
In [ ]:
In [ ]:
In [1246]:
Gaussian_length3 = 10000
Gaussian_LIBOR = pd.DataFrame()
LIBOR_top_limit = 20
LIBOR_bottom_limit = 0.1
number_of_simulated_days = 2**13+1
start_time = datetime.now().strftime("%I:%M%p")
for i in range(0,Gaussian_length3):
print("Having started at "+str(start_time))
print("i="+str(i))
Gaussian_LIBOR[i] = gaussian_attempt(LIBOR_top_limit, LIBOR_bottom_limit, df.LN_LIBOR.mean(), df.LN_LIBOR.std(), df.Price_LIBOR[0], number_of_simulated_days)
clear_output()
print("Having started at "+str(start_time))
Gaussian_LIBOR.to_csv('Gaussian_LIBOR.csv')
In [1247]:
plt.figure(figsize=(24,10))
for i in range(0,Gaussian_length3):
K=13
plt.title(str(Gaussian_length3) + " Gaussian generated price paths")
plt.xticks(np.arange(0, 2**(K)+1, 2**(K-3)))
plt.xlabel("Time\n(days)")
plt.ylabel('Price level')
plt.grid(b=True)
plt.plot(Gaussian_LIBOR[i], color="cornflowerblue", linewidth=0.1, alpha=0.2)
plt.show()
In [1248]:
# Here we calculate the LN differences
Gaussian_diff_attempt3 = [np.log(Gaussian_LIBOR[x]).diff() for x in range(0,Gaussian_length3)]
Gaussian_diff_attempt3 = pd.DataFrame(Gaussian_diff_attempt3).transpose()
Gaussian_diff_attempt3 = Gaussian_diff_attempt3[:-1]
Gaussian_diff_attempt3
Gaussian_diff_attempt3.to_csv('Gaussian_diff_LIBOR.csv')
In [71]:
for i in range(0,5):
plt.figure(figsize=(24,5))
plt.title("Price increments for simulation " + str(i+1))
plt.xticks(np.arange(0, 2**(13)+1, 2**(13-3)))
plt.xlabel("Time\n(days)")
plt.grid(b=True)
plt.ylabel('Change\n(%)')
axes = plt.gca()
axes.set_ylim([-0.25,0.25])
plt.plot(Gaussian_diff_attempt3[i], color="forestgreen", linewidth=0.5)
plt.show()
In [95]:
plt.grid(b=True)
plt.title("Histogram of Gaussian simulated final prices\nafter 8192 days for LIBOR")
plt.hist(Gaussian_LIBOR.iloc[-1], bins=40, color="cornflowerblue")
plt.xlabel("Final yield (%)")
plt.ylabel("Frequency")
plt.axvline(Gaussian_LIBOR.iloc[-1].quantile(0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(Gaussian_LIBOR.iloc[-1].quantile(0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(Gaussian_LIBOR.iloc[-1].quantile(0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(Gaussian_LIBOR.iloc[-1].quantile(0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(Gaussian_LIBOR.iloc[-1].quantile(0.95), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_LIBOR.iloc[-1], linestyle='--', color="black", label="real final price", linewidth=0.7)
plt.legend()
plt.show()
In [97]:
plt.grid(b=True)
plt.title("Mean yield change for \n" + str(Gaussian_length3) + " Gaussian generated price paths\nfor LIBOR")
plt.hist(Gaussian_diff_attempt3.mean(), bins=30, color='cornflowerblue')
plt.axvline(df.LN_LIBOR.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_LIBOR.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(Gaussian_length3))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((Gaussian_diff_attempt3.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of yield changes for \n" + str(Gaussian_length3) + " Gaussian generated price paths\nfor LIBOR")
plt.hist(Gaussian_diff_attempt3.std(), bins=30, color='cornflowerblue')
plt.axvline(df.LN_LIBOR.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_LIBOR.std(), 5)))
plt.ylabel("Frequency per " + str(Gaussian_length3))
plt.xlabel("Standard Deviation")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our 10,000 simulations.\n")
(Gaussian_diff_attempt3.std()).describe()
Out[97]:
In [181]:
print("The mean is obviously very variable - we can simulate both downward and upward markets.\n\nWe can see that the standard deviation is also quite variable, but it is mostly near to what we expect.")
In [182]:
print("Here's the kurtosis for the real data for Norway.")
stats.kurtosis(df.LN_NORWAY[1:-1])
Out[182]:
In [82]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length3) + " Gaussian generated price paths\nfor Norway")
plt.hist(stats.kurtosis(Gaussian_diff_attempt3), bins=100, color='cornflowerblue')
plt.axvline(stats.kurtosis(df.LN_LIBOR[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_LIBOR[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt3 = stats.kurtosis(Gaussian_diff_attempt3)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length3))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [92]:
plt.figure(figsize=(16,4))
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length3) + " Gaussian generated price paths\nfor Norway")
plt.hist(stats.kurtosis(Gaussian_diff_attempt3), bins=100, color='cornflowerblue')
# plt.axvline(stats.kurtosis(df.LN_LIBOR[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_LIBOR[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt3 = stats.kurtosis(Gaussian_diff_attempt3)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length3))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [102]:
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length3) + " Gaussian generated price paths\nfor Norway")
plt.hist(stats.kurtosis(Gaussian_diff_attempt3), bins=100, color='cornflowerblue')
# plt.axvline(stats.kurtosis(df.LN_LIBOR[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_LIBOR[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
kurtosis_array_Gaussian_diff_attempt3 = stats.kurtosis(Gaussian_diff_attempt3)
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3, 0.05), linestyle='--', color="purple", label="5th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.25), linestyle='--', color="blue", label="25th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.5), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.75), linestyle='--', color="gold", label="75th percentile")
plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.95), linestyle='--', color="red", label="95th percentile")
plt.ylabel("Frequency per " + str(Gaussian_length3))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [83]:
# Here we test how well the simulations fit the original data, using the Kolmogorov-Smirnov test to test the log price changes
ks_test_Gaussian_LIBOR = [0 for x in range(0, Gaussian_length3)]
for i in range(0,Gaussian_length3):
ks_test_Gaussian_LIBOR[i] = stats.ks_2samp(Gaussian_diff_attempt3[i], df.LN_LIBOR).pvalue
plt.grid(b=True)
plt.hist(ks_test_Gaussian_LIBOR, bins=20, color='cornflowerblue')
plt.title("Histogram of Kolmogorov-Smirnov P-values\nfor " + str(Gaussian_length3) + " Gaussian simulations for Norway ")
plt.axvline(0.05, linestyle='--', color="black", label="5% significance")
plt.xlabel("P-Value")
plt.ylabel("Frequency per " + str(Gaussian_length3))
plt.legend()
plt.show()
plt.grid(b=True)
plt.plot(sorted(ks_test_Gaussian_LIBOR), color='cornflowerblue')
plt.title("Sorted Kolmogorov-Smirnov P-values\nfor " + str(Gaussian_length3) + " Gaussian simulations for Norway ")
plt.axhline(0.05, linestyle='--', color="black", label="5% significance")
plt.axvline(sum(x<0.05 for x in ks_test_Gaussian_LIBOR), linestyle='--', color="forestgreen", label="Total rejections = "+str(sum(x<0.05 for x in ks_test_Gaussian_LIBOR)))
plt.xlabel("Index ouf of " + str(Gaussian_length3))
plt.ylabel("P-Value")
plt.legend()
plt.show()
print("Count of P<0.05 = "+ str(sum(x<0.05 for x in ks_test_Gaussian_LIBOR)) + " failures out of " + str(Gaussian_length3) + ".")
In [ ]:
In [103]:
print("Let's get all the kurtosis figures onto one graph")
In [107]:
# NORWAY GRAPHS
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length) + " Gaussian generated price paths\nfor Norway")
plt.hist(stats.kurtosis(Gaussian_diff_attempt), bins=100, color='orangered', alpha = 0.5, label = "Norway")
# plt.axvline(stats.kurtosis(df.LN_NORWAY[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_NORWAY[1:-1]), 3)))
plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
# kurtosis_array_Gaussian_diff_attempt = stats.kurtosis(Gaussian_diff_attempt)
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt, 0.05), linestyle='--', color="purple", label="5th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.25), linestyle='--', color="blue", label="25th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.5), linestyle='--', color="forestgreen", label="50th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.75), linestyle='--', color="gold", label="75th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt,0.95), linestyle='--', color="red", label="95th percentile")
# SWEDEN GRAPHS
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length2) + " Gaussian generated price paths\nfor Sweden")
plt.hist(stats.kurtosis(Gaussian_diff_attempt2), bins=100, color='goldenrod', alpha = 0.5, label = "Sweden")
# plt.axvline(stats.kurtosis(df.LN_SWEDEN[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_SWEDEN[1:-1]), 3)))
# plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
# kurtosis_array_Gaussian_diff_attempt2 = stats.kurtosis(Gaussian_diff_attempt2)
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2, 0.05), linestyle='--', color="purple", label="5th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.25), linestyle='--', color="blue", label="25th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.5), linestyle='--', color="forestgreen", label="50th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.75), linestyle='--', color="gold", label="75th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt2,0.95), linestyle='--', color="red", label="95th percentile")
# LIBOR GRAPHS
plt.grid(b=True)
plt.title("Kurtosis of price changes for \n" + str(Gaussian_length3) + " Gaussian generated price paths\nfor Norway")
plt.hist(stats.kurtosis(Gaussian_diff_attempt3), bins=100, color='cornflowerblue', alpha = 0.5, label = "LIBOR")
# plt.axvline(stats.kurtosis(df.LN_LIBOR[1:-1]), linestyle='--', color="black", label="real kurtosis\n= "+str(round(stats.kurtosis(df.LN_LIBOR[1:-1]), 3)))
# plt.axvline(0, linestyle='--', color="aqua", label="Gaussian kurtosis = 0")
# kurtosis_array_Gaussian_diff_attempt3 = stats.kurtosis(Gaussian_diff_attempt3)
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3, 0.05), linestyle='--', color="purple", label="5th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.25), linestyle='--', color="blue", label="25th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.5), linestyle='--', color="forestgreen", label="50th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.75), linestyle='--', color="gold", label="75th percentile")
# plt.axvline(np.quantile(kurtosis_array_Gaussian_diff_attempt3,0.95), linestyle='--', color="red", label="95th percentile")
# PLT.SHOW
plt.ylabel("Frequency per " + str(Gaussian_length3))
plt.xlabel("Kurtosis")
plt.legend()
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [180]:
print("Our graphs before were a bit messy.\n\nLet's make them messier by adding a direct comparison to the Gaussian on the same graph.")
In [179]:
plt.grid(b=True)
plt.title("Histogram of simulated final prices\nafter 8192 days for Norway")
plt.hist(Gaussian_Norway.iloc[-1], bins=40, color="aqua", alpha = 0.2, label = "Gaussian comparison")
plt.hist(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191]), bins=40, color="orangered")
plt.xlabel("Final price (USD/NOK)")
plt.ylabel("Frequency")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.05)), linestyle='--', color="purple", label="5th percentile")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.25)), linestyle='--', color="blue", label="25th percentile")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.5)), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.75)), linestyle='--', color="gold", label="75th percentile")
plt.axvline(df.Price_NORWAY[0]*np.exp(attempt.iloc[8191].quantile(0.95)), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_NORWAY.iloc[-1], linestyle='--', color="black", label="real final price", linewidth=0.7)
plt.legend()
plt.show()
In [202]:
plt.grid(b=True)
plt.title("Mean price change for \n" + str(attempt_length) + " MMAR generated price paths\nfor NORWAY")
plt.hist(diff_attempt.mean(), bins=30, color='orangered')
plt.hist(Gaussian_diff_attempt.mean(), bins=60, color='aqua', alpha = 0.2, label = "Gaussian comparison")
plt.axvline(df.LN_NORWAY.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_NORWAY.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(attempt_length))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((diff_attempt.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of price changes for \n" + str(attempt_length) + " MMAR generated price paths\nfor NORWAY")
plt.hist(diff_attempt.std(), bins=30, color='orangered')
plt.axvline(df.LN_NORWAY.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_NORWAY.std(), 5)))
plt.ylabel("Frequency per " + str(attempt_length))
plt.xlabel("Standard Deviation")
plt.hist(Gaussian_diff_attempt.std(), bins=15, color='aqua', alpha = 0.2, label = "Gaussian comparison")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our 10,000 simulations.\n")
(diff_attempt.std()).describe()
Out[202]:
In [161]:
plt.grid(b=True)
plt.title("Histogram of simulated final prices\nafter 8192 days for Sweden")
plt.hist(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191]), bins=40, color="goldenrod")
plt.hist(Gaussian_Sweden.iloc[-1], bins=40, color="aqua", alpha = 0.2, label = "Gaussian comparison")
plt.xlabel("Final price (SEK)")
plt.ylabel("Frequency")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.05)), linestyle='--', color="purple", label="5th percentile")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.25)), linestyle='--', color="blue", label="25th percentile")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.5)), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.75)), linestyle='--', color="gold", label="75th percentile")
plt.axvline(df.Price_SWEDEN[0]*np.exp(attempt2.iloc[8191].quantile(0.95)), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_SWEDEN.iloc[-1], linestyle='--', color="black", label="real final price")
plt.legend()
plt.show()
In [204]:
plt.grid(b=True)
plt.title("Mean price change for \n" + str(attempt2_length) + " MMAR generated price paths\nfor SWEDEN")
plt.hist(diff_attempt2.mean(), bins=30, color='goldenrod')
plt.hist(Gaussian_diff_attempt2.mean(), bins=30, color='aqua', alpha = 0.2, label = "Gaussian comparison")
plt.axvline(df.LN_SWEDEN.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_SWEDEN.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(attempt2_length))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((diff_attempt2.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of price changes for \n" + str(attempt2_length) + " MMAR generated price paths\nfor SWEDEN")
plt.hist(diff_attempt2.std(), bins=30, color='goldenrod')
plt.hist(Gaussian_diff_attempt2.std(), bins=30, color='aqua', alpha = 0.2, label = "Gaussian comparison")
plt.axvline(df.LN_SWEDEN.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_SWEDEN.std(), 4)))
plt.ylabel("Frequency per " + str(attempt2_length))
plt.xlabel("Standard Deviation")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our " + str(attempt2_length) + " simulations.\n")
(diff_attempt2.std()).describe()
Out[204]:
In [169]:
plt.grid(b=True)
plt.title("Histogram of simulated final yields\nafter 8192 days for LIBOR")
plt.hist(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191]), bins=40, color="cornflowerblue")
plt.hist(Gaussian_LIBOR.iloc[-1], bins=40, color="aqua", alpha = 0.2, label = "Gaussian comparison")
plt.xlabel("Final yield (%)")
plt.ylabel("Frequency")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.05)), linestyle='--', color="purple", label="5th percentile")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.25)), linestyle='--', color="blue", label="25th percentile")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.5)), linestyle='--', color="forestgreen", label="50th percentile")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.75)), linestyle='--', color="gold", label="75th percentile")
plt.axvline(df.Price_LIBOR[0]*np.exp(attempt3.iloc[8191].quantile(0.95)), linestyle='--', color="red", label="95th percentile")
plt.axvline(df.Price_LIBOR.iloc[-1], linestyle='--', color="black", label="real final price")
plt.legend()
plt.show()
In [172]:
plt.grid(b=True)
plt.title("Mean yield change for \n" + str(attempt3_length) + " MMAR generated price paths\nfor 12-month LIBOR")
plt.hist(diff_attempt3.mean(), bins=30, color='cornflowerblue')
plt.hist(Gaussian_diff_attempt3.mean(), bins=30, color='aqua', alpha = 0.2, label = "Gaussian comparison")
plt.axvline(df.LN_LIBOR.mean(), linestyle='--', color="black", label="real mean\n= "+str(round(df.LN_LIBOR.mean(), 6)))
plt.axvline(0, linestyle='--', color="forestgreen", label="mean = zero")
plt.ylabel("Frequency per " + str(attempt3_length))
plt.xlabel("Mean")
plt.legend()
plt.show()
print((diff_attempt3.mean()).describe())
plt.grid(b=True)
plt.title("Standard deviation of yield changes for \n" + str(attempt3_length) + " MMAR generated price paths\nfor 12-month LIBOR")
plt.hist(diff_attempt3.std(), bins=30, color='cornflowerblue')
plt.hist(Gaussian_diff_attempt3.std(), bins=30, color='aqua', alpha = 0.2, label = "Gaussian comparison")
plt.axvline(df.LN_LIBOR.std(), linestyle='--', color="black", label="real standard deviation\n= "+str(round(df.LN_LIBOR.std(), 5)))
plt.ylabel("Frequency per " + str(attempt3_length))
plt.xlabel("Standard Deviation")
plt.legend()
plt.show()
print("\n\nBelow is a description of the standard deviations of our " + str(attempt3_length) + " simulations.\n")
print((diff_attempt3.std()).describe())
In [ ]:
In [ ]: